
#include <msOscillatorsAssembly.h>


using namespace atomism;
boost::shared_ptr<msOscillator> castInOscill( boost::shared_ptr<msTreeMapper> ptr ){ return boost::static_pointer_cast<msOscillator>(ptr);  
										    }

int main(int argc, char **argv) {
 
 boost::shared_ptr<msUnitsManager> unitMng = msUnitsManager::New();
 
 unitMng->set( "Angstrom Degree amu kcal/mol ps" );
 
 boost::shared_ptr<msOscillatorsAssembly> oscillsHO = msOscillatorsAssembly::New(unitMng);
 
 oscillsHO->addOscillator( castInOscill( msHarmonicOscillator::New("cm-1")->set(1689.00) ) );
 oscillsHO->addOscillator( castInOscill( msHarmonicOscillator::New("cm-1")->set(3828.44) ) );
 oscillsHO->addOscillator( castInOscill( msHarmonicOscillator::New("cm-1")->set(3934.84) ) );
 oscillsHO->getParameters()->setParameter("Emax",100./3, *unitMng);
 oscillsHO->computeDOS();

 double Qho = oscillsHO->Q(1500);
 double QhoDos = oscillsHO->QfromDOS(1500);
 
 boost::shared_ptr<msOscillatorsAssembly> oscillsAO = msOscillatorsAssembly::New(unitMng);

 oscillsAO->addOscillator( castInOscill( msAnharmonicOscillatorP2::New("cm-1")->set(1689.00,"-15.88 -22.18 -29.2")   ));
 oscillsAO->addOscillator( castInOscill( msAnharmonicOscillatorP2::New("cm-1")->set(3828.44,"-22.18 -46.18 -172.14") ));
 oscillsAO->addOscillator( castInOscill( msAnharmonicOscillatorP2::New("cm-1")->set(3934.84,"-29.2 -172.14 -60.09")  ));
 oscillsAO->getParameters()->setParameter("Emax",100./3, *unitMng);
 oscillsAO->computeDOS();

 double Qao = oscillsAO->Q(1500);
 double QaoDos = oscillsAO->QfromDOS(1500);
 
 fprintf(stdout,"Qho(1500)=%g ; ",    Qho);
 fprintf(stdout,"QhoDos(1500)=%g ; ", QhoDos);
 fprintf(stdout,"Qao(1500)=%g ; ",    Qao);
 fprintf(stdout,"QaoDos(1500)=%g ; ", QaoDos);
 fprintf(stdout,"\n");
 return 0;
 }


